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ABSTRACT 

Motivated by the nonlinear star formation efficiency found in recent numerical 
simulations by a number of workers, we perform high-resolution adaptive mesh refine¬ 
ment simulations of star formation in self-gravitating turbulently driven gas. As we 
follow the collapse of this gas, we find that the character of the flow changes at two 
radii, the disk radius rg, and the radius r* where the enclosed gas mass exceeds the 
stellar mass. Accretion starts at large scales and works inwards. In line with recent 
analytical work, we find that the density evolves to a fixed attractor, p{r,t) —i p(r), 
for ry < r < mass flows through this structure onto a sporadically gravitationally 
unstable disk, and from thence onto the star. In the bulk of the simulation box we 
find that the random motions vt ~ with p ~ 0.5, in agreement with Larson’s size- 
linewidth relation. In the vicinity of massive star forming regions we find p ~ 0.2 — 0.3, 
as seen in observations. For r < r*, vp increases inward, with p = —1/2. Finally, we 
find that the total stellar mass M* (t) ~ in line with previous numerical and analytic 
work that suggests a nonlinear rate of star formation. 

Key words: galaxies: star clusters: general - galaxies: star formation - stars: forma¬ 
tion - turbulence 


1 INTRODUCTION 

The star formation time on galactic scales is long when com¬ 
pared to the dynamical time. Kennicutt (1998) expressed 
this in the form 

S, = eSgTQYN (1) 

where S* is the star formation rate per unit area, Eg is 
the gas surface density, tdyn is the local dynamical time, 
and e = 0.017 is the efficiency factor. Naively, if the gas 
self-gravity dominates the dynamics, e ~ 1, so the low effi¬ 
ciency of star formation is surprising. More recent work has 
refined this and similar relations in regard to its dependence 
on molecular gas (Bigiel et al. 2008) and by taking into ac¬ 
count the error distributions of both S* and Eg (Shetty et al. 
2013), but the best current estimates of the efficiency of star 
formation on galactic scales remains low. 

Whether this low efficiency applies to scales comparable 
to giant molecular clouds, with radii of order 100 pc, is de¬ 
bated in the literature. Heiderman et al. (2010), Lada et al. 
(2010), Wu et al. (2010), and Murray (2011) find efficien¬ 
cies a factor of ten or more larger, while Krumholz & Tan 
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(2007) and Krumholz et al. (2012a) find e ~ 0.01. On these 
small scales, observations also suggest that the efficiency is 
not universal, but instead varies over two to three orders of 
magnitude (e.g. Mooney & Solomon 1988; Lee et al. 2016). 

There are a number of explanations for the low star for¬ 
mation rate, on either small or large scales (although they 
may not be necessary for the former!). These include turbu¬ 
lent pressure support (Myers & Fuller 1992), support from 
magnetic fields (Strittmatter 1966; Mouschovias 1976), and 
stellar feedback (e.g. Dekel & Silk 1986). Numerical experi¬ 
ments investigating the first two effects suggest that neither 
turbulence nor magnetic support is sufficient to reduce the 
rate of star formation to e « 0.02 on small scales (Wang 
et al. 2010; Cho & Kim 2011; Padoan & Nordlund 2011; 
Krumholz et al. 2012b; Myers et al. 2014). Feedback from 
radiative effects and protostellar jets and winds may be able 
to explain the low star formation rate, but the impact of 
these forms of stellar feedback remains uncertain despite re¬ 
cent progress (Wang et al. 2010; Myers et al. 2014; Federrath 
2015). 

Until very recently, galaxy-scale or larger (cosmologi¬ 
cal) simulations were not able to reproduce the Kennicutt- 
Schmidt relation. Nor did the cosmological runs reproduce 
correctly the mass of stars in galaxies of a given halo mass, 
despite including supernova and other forms of feedback. 
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e.g., Guo et al. (2010); Governato et al. (2010); Piontek & 
Steinmetz (2011). To overcome this low resolution driven 
problem, Hopkins et al. (2011, 2012) performed high reso¬ 
lution (few parsec spatial, few hundred solar mass particle 
masses) simulations of isolated galaxies, modeling both ra¬ 
diative and supernovae feedback (among other forms). They 
recovered the Kennicutt-Schmidt relation, a result that they 
showed was independent of the small-scale star formation 
law that they employed. The simulations in the second pa¬ 
per also generated galaxy scale outflows or winds, removing 
gas from the disk, thus making it unavailable for star forma¬ 
tion. When the feedback was turned off, the star formation 
rate soared, demonstrating that in the simulations at least, 
feedback was crucial to explaining the Kennicutt-Schmidt 
relation, and the outflows. Simulations including supernovae 
but lacking the radiative component of the feedback did not 
exhibit strong winds and so overproduced stars. 

Gosmological simulations employing unresolved (or 
“sub-grid”) models for both radiative and supernovae feed¬ 
back are now able to reproduce the halo-mass/stellar mass 
relation (e.g., Aumer et al. 2013; Hopkins et al. 2014; Agertz 
& Kravtsov 2015). Again, these simulations require stellar 
feedback to drive the winds that remove gas from the disk, 
so as to leave the observed mass of stars behind. 

Lee et al. (2015) emphasized that the star formation ef¬ 
ficiency on parsec scales is nonlinear in time, i.e., e oc t —> 
M* (X on small scales, where M* is the total stellar mass. 
Magnetic fields slowed the initial star formation rate some¬ 
what, but did not change the M*(t) ~ scaling. Using a 
detailed numerical simulation, they showed that this non¬ 
linear star formation rate is driven by the properties of col¬ 
lapsing regions. In particular, they demonstrated that the 
turbulent velocity near or in collapsing regions follows dif¬ 
ferent scaling relations than does turbulence in the global 
environment, which follows Larson’s law, vtO ~ (Lar¬ 
son 1981). They also showed that the density PDF is not 
log-normal, but rather develops a power law to high density. 
This latter result was hinted at by Klessen (2000) and shown 
convincingly, as well as explained, by Kritsuk et al. (2011). 

The increasing rate of star formation found by Lee et al. 
(2015) is important in that it may provide an explanation for 
the observed range in star formation rates on small scales. It 
suggests that the star formation rates on small scales vary in 
part because of the age of the star forming region; slow star 
forming regions, with very low instantaneous efficiencies, will 
ramp up their stellar production over time. If this result can 
be firmly established, it will highlight the need for a form of 
very rapid feedback. In particular, since the dynamical time 
in massive star forming regions is much smaller than the 
time delay of ~ 4Myrs between the start of star formation 
and the first supernovae, rapid star formation on small scales 
would have to be halted by some form of feedback other than 
supernovae. 

The simulations of Lee et al. (2015) explicate the link 
between the rate of star formation with the gravitational 
collapse of high density regions, which is an analytically well 
studied problem. An early model of Shu (1977) estimated the 
accretion rate onto stars by assuming that stars form from 
hydrostatic cores supported by thermal gas pressure. The 
accretion rate in his model was independent of time, given 
by M = moc?s/G, where Cs = {kbT/ii)^G jg the sound speed 
in molecular gas, and mo = 0.975. Shu (1977) predicted a 


maximum accretion rate of ~ 2 x 1O“®M0 yr“^, which is too 
small to explain the origin of massive (M* ~ 50 — IOOMq) 
O stars, which have lifetimes < 4 x 10® yrs 

Myers & Fuller (1992) overcame the difficulty with slow 
accretion rates by adopting the turbulent speed in lieu of 
the sound speed (see also McLaughlin & Pudritz (1997) and 
McKee & Tan (2003)). In doing so they were able to replace 
the slower signal speed of sound with the faster turbulent 
speed. However, they continued to assume the initial con¬ 
dition was that of a hydrostatic core that is supported by 
turbulent pressure. They also assumed that the turbulence 
is static and unaffected by the collapse. 

Gollectively, these models, (Shu 1977; Myers & Fuller 
1992; McLaughlin & Pudritz 1997; McKee & Tan 2003), are 
referred to as inside-out collapse models; the collapse starts 
at small radii (formally at r = 0 in the analytic models) and 
works its way outward, at the assumed propagation speed 
(Cs or vtO). At any given time, the infall velocity and mass 
accretion rate both decrease with increasing radius r. The 
analytic models assume the existence of a self-similarity vari¬ 
able X = r/vt, where u = Cs in Shu (1977) or the turbulent 
velocity vtO in Myers & Fuller (1992); McLaughlin & Pu¬ 
dritz (1997); McKee & Tan (2003). These models predict ve¬ 
locity and mass accretion profiles very different than those 
seen in the simulations of Lee et al. (2015). 

Motivated by this discrepancy, Murray & Chang (2015), 
hereafter MC15, developed a 1-D model of spherical col¬ 
lapse that treats the turbulent velocity, vt, as a dynamical 
variable and does not assume that the initial condition is a 
hydrostatically supported region. They used the results of 
Robertson & Goldreich (2012) on compressible turbulence; 
the evolution of the turbulent velocity in a collapsing (or ex¬ 
panding) region is described well by the following equation: 


dVT 


Ur 


dVT 

dr 


+ 



VtUt 

r 


= 0 


( 2 ) 


The first two terms are the Lagrangian derivative, and Ur 
is the radial infall velocity. The first term in the brackets 
describes the turbulent driving produced by the infall, while 
the second is the standard expression for the turbulent decay 
rate; is a dimensionless constant of order unity. 

MG15 used this in place of an energy equation. To¬ 
gether with the equations for mass continuity and momen¬ 
tum, equation (2) gives a closed set of equations that can be 
solved in spherical symmetry numerically. In addition, they 
were able to analytically show that the results of their cal¬ 
culations gave density and velocity profiles that appear to 
be in line with both recent numerical calculations (Lee et al. 
2015) and observations (e.g., Caselli & Myers 1995; Plume 
et al. 1997). 

To summarize, MC15’s major results were: 


• The gravity of the newly formed star introduces a phys¬ 
ical scale into the problem, which MC15 called the stellar 
sphere of influence, r*. This is an idea familiar from galactic 
dynamics. The radius r* is where the local dynamics transi¬ 
tions from being dominated by the mass of the gas to being 
dominated by the mass of the star. As a result,the charac¬ 
ter of the solution, in particular that of the velocity, differs 
dramatically between r < r* and r > r*. The existence of 
this physical scale modifies the form of the self-similarity on 
which inside-out theories rely. 
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• The small scale density profile is an attractor solution. 
MC15 showed numerically and argued analytically that at 
small scales, the density profile is an attractor solution. In 
particular, MC15 showed the density profile asymptotes to: 


p(r, t) = { 



p{ro,t) 



hp 


r < rt,{t) 
1.6 — 1.8 r > r*(t). 


(3) 


where ro is some fiducial radius. 

• The existence of r, implies that the infall and turbulent 
velocities have different scaling for r < r* and r > r*. In 
particular, MC15 showed 


Ur{r,t),VT{r,t) oc < 


GM4t) 

r 


GM{r,t) ^ ^ 0.2 


r < rt,{t) 
r > rt{t), 


(4) 


Thus the scaling of the turbulent velocity differs from that 
predicted by Larson’s law (oc inside the sphere of influ¬ 
ence. In other words, the turbulent velocity in massive star 
forming regions will deviate from Larson’s law, which has 
long been observed, but without theoretical explanation. 

• The stellar mass increases quadratically with time. This 
result arises naturally from the attractor solution nature of 
the density profile at small r. Equation (3), and the scaling 
with Keplerian velocity for the turbulent and infall velocities 
at small r, Equation (4). 

The mass accretion rate: 


M(r, t) 


‘iirB?p{R)ur{r,t),~ tC r < r* 

4:TtR^ p{R)ur{r,t) ~ T r > r*. 


(5) 


MClS’s predictions for r < r* could not be checked us¬ 
ing the simulations of Lee et al. (2015) as those fixed grid 
simulations were too coarse. In this paper, we study the col¬ 
lapse of gas and formation of stars in a turbulent GMC using 
roughly a dozen high resolution adaptive mesh refinement 
(AMR) simulations in this paper. 

We employ large-scale (16 pc) hydrodynamic AMR sim¬ 
ulations of star-forming clouds with continuously driven su¬ 
personic turbulence. The initial conditions for our simula¬ 
tions are exactly the same as the FLASH simulations in Lee 
et al. (2015). 

If the equations are non-dimensionalized, two dimen¬ 
sionless variables appear, the Mach number A4 and the virial 
parameter Ovir = {5/3)vf/GMbox, e.g., Mihalas & Mihalas 
(1984). We want to model massive star forming regions in 
the Milky Way, so we choose the Mach Number A4 = 9 
and the virial parameter avir = 1.9 respectively. In addi¬ 
tion, we choose the size of the box L = 16pc, and the sound 
speed Cs = 0.264 km s“^, so that the turbulent velocity lies 
approximately on the observed size-line width relation, Lar¬ 
son’s Law. These choices fix both the density and the mass 
scale. 

The simulations described in this paper disregard sev¬ 
eral physical effects. We do not include radiative, stel¬ 
lar wind, or proto-stellar jet feedback. While the feedback 
physics we neglect can have significant effects on both the 


rate of star formation and the initial mass function (IMF), 
we aim to address the role the random motions captured 
by the Reynolds stress play in the dynamics of gravitational 
collapse in turbulent fluids. 

Our equation of state is that of an isothermal gas. It is 
possible, and even likely, that thermal effects play a role in 
setting the initial mass function of stars, e.g Larson (2005). 
With this in mind, we relegate the discussion of the IMF to 
an appendix, as the details are unlikely to be reliable. 

This paper is organized as follows. In Section 2 we de¬ 
scribe our numerical methods and simulation setup. In Sec¬ 
tion 3 we present and analyze the results of our simulations. 
In particular, we make detailed comparisons with the re¬ 
sults of MC15. We discuss our results and compare them to 
previous work in Section 4 


2 DETAILED SIMULATIONS OF TURBULENT 

COLLAPSE 

Most of the simulations described here use the adaptive 
mesh refinement code FLASH ver. 4.0.1 (Fryxell et al. 2000; 
Dubey et al. 2008) to model self-gravitating, hydrodynamic 
turbulence in isothermal gas with three-dimensional (3D) 
periodic grids and a minimum of 8 levels of refinement on a 
root grid of 128^, giving an effective resolution of 32,768®. 
Following Lee et al. (2015) our FLASH runs use the Harten- 
Lax-van Leer-Contact Riemann solver and an unsplit solver 
(Lee et al. 2009). We have also used the RAMSES code 
(Teyssier 2002), but unless explicitly stated otherwise, the 
results below come from FLASH simulations. 

As just mentioned, we start with a box with the physical 
length set to L = 16 pc using periodic boundary conditions. 
The initial mass density is p = 3 x 10“®®gcm“® (number 
density n ~ 100cm“®), corresponding to a mean free-fall 
time Tff « 3.8Myrs; the total mass in the box is M « 
18,000M©. The sound speed is set to Ca = 0.264kms“®. 
We use pure molecular hydrogen in this simulation so the 
ambient temperature T « 17K. 

To initialize our simulations, we drive turbulence by 
applying a large scale (1 < kL < 2, corresponding to 
1.3 — 2.7 pc) fixed solenoidal acceleration field as a momen¬ 
tum source term. We use solenoidal driving because it is 
known that compressive turbulence increases the star for¬ 
mation rate compared to solenoidal driving (Federrath et al. 
2008). We apply this field in the absence of gravity and 
star particle formation for 3 dynamical times until a sta¬ 
tistical steady state is reached. The resulting Mach number 
is Af = 9, i.e a turbulent velocity of vt = 2.37kms“®. 

Stirring the initial turbulence using a fixed driving field 
is a technique used by a number of workers in the field 
(Padoan & Nordlund 2011; Collins et al. 2011). Other groups 
initialize the turbulence by initializing the velocity field with 
Caussian random perturbations having some assumed power 
spectrum (Myers et al. 2014; Skinner & Ostriker 2015). 
While neither of the resulting velocity fields are generated 
the way the turbulence in the interstellar medium (ISM) of 
our Calaxy is, the stirring allows one to perform simulations 
which have nontrivial initial density structures and velocity 
fields that are at least reminiscent of those inferred from 
observations of the interstellar medium of our Calaxy. 

Federrath & Klessen (2012) use a time varying driving 
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field to produce random motions. They argue that a time- 
varying driving field allows one to avoid large spatial cor¬ 
relations that would result from a fixed driving field acting 
for a time longer than the dynamical time of the simulation 
box. In our simulations we do not run for longer than a box 
dynamical time after turning on star formation. We run for 
600,000yrs, about 0.16 dynamical times, after the first star 
forms. The limiting factor on the length of the runs was our 
available compute time. Hence, we do not expect the large 
scale turbulent flow to vary much over such a short time. 
In addition, there is some evidence (Federrath et al. 2010a) 
that the results of turbulent driving are not sensitive to the 
exact large-scale mechanism. 

This fully developed turbulent state is the initial condi¬ 
tion to which we add self-gravity and star particle formation 
for our star formation experiments. We enable AMR to fol¬ 
low the collapse of overdense regions. Even after turning on 
star formation, we continue to drive the large scale fixed 
solenoidal acceleration field. 

To follow these collapsing regions, we have implemented 
an algorithm for mesh refinement in these simulations, sim¬ 
ilar to that of Federrath et al. (2010b). In supersonically 
turbulent flows, certain regions rapidly increase in density. 
For a given density and temperature, or sound speed, regions 
larger than the Jeans length 



are prone to gravitational collapse. Our base grid’s resolu¬ 
tion of N^oot = 128® gives a cell length of 1.25xl0“^pc which 
is sufficient to resolve the Jeans length for the mean density. 

In most of our simulations, the AMR grid is refined 
when the Truelove et al. (1997) criterion 

\j < NjAx, (7) 


is met. In this expression Aa; is the cell length, and Nj is an 
integer; Truelove et al. (1997) found that in order to avoid 
artificial fragmentation, one requires Nj > 4. 

This corresponds to a condition on the density 


Po 


45 • 4' 


h Nioot 

V 128 



0.265 km s 


2 ^3 X 10 ®®gcm 
po 


( 8 ) 


where I is the refinement level, with I — 0 corresponding to 
the root grid. When this density condition is met the local 
grid is refined by a factor of 2, provided that the maximum 
refinement level has not been reached. When the transisition 
to the maximum refinement level is triggered i.e. when I goes 
from 7 to 8 (the maximum refinement level), the density 
contrast is p/po ~ 10®. 

In the Appendix we describe a number of test simula¬ 
tions in which we refined the grid when Nj — 4, 8,16 or 32 
(Federrath et al. 2011). We show that many of the quanti¬ 
ties in our runs, including the density and the mass accretion 
rates, are converged for Nj — 4. 

The maximum dynamic range is a little larger than 6 or¬ 
ders of magnitude, because we allow the density to increase 
further before forming star particles. When the Truelove cri¬ 
terion is exceeded by a factor of three at the highest refine¬ 
ment level, the excess mass in a cell is transferred either to 


a newly created star particle or to a star particle whose ac¬ 
cretion radius includes the cell. The factor of three allows 
only the highest density regions to form star particles. It is 
inspired by the work of Padoan & Nordlund (2011) whose 
sink particle formation criteria of 8000 x mean density is a 
factor of 3-4 above the Truelove criteria at their highest reso¬ 
lution of 1000®. Additionally, the 3 cells immediately around 
a star particle can rise above this density criterion. This is 
done so that we do not form star particles within 2 cells of 
each other. Instead these close surrounding cells can only 
accrete onto the previously formed star particle. We should 
also note that like our previous work in Lee et al. (2015), our 
star particle creation prescription is different from the pre¬ 
scription of Federrath et al. (2010b) where additional checks 
are performed; in the appendix we present the results of 
runs in which we used these additional checks, finding that 
they do not affect the scaling of the stellar mass, or the 
dynamics of the infall. 

To calculate the gravity, we use the same algorithm as 
described in Lee et al. (2015), which we now briefly describe. 
To compute the self gravity on gas, we first map star parti¬ 
cles to the grid and then use a multi-grid Poisson solver (see 
Ricker 2008), coupled with a fast-Fourier transform (FFT) 
solution on the root grid, to solve for gravity. To compute 
the gravitational acceleration on the star particles, we first 
compute the particle-particle forces using a direct N-body 
calculation. To compute the particle-gas forces, we use the 
same multigrid solver (with root grid FFT) on the grid, but 
with the star particle unmapped. As a result, two large scale 
gravity solutions (one with and one without mapped star 
particles) must be found per timestep as opposed to one. 
This allows us to avoid the computationally expensive task 
of computing gas-star particle forces via direct summation. 
As discussed in Lee et al. (2015), this splitting of particle- 
particle and particle-gas/gas-particle forces does not strictly 
obey Newton’s second law, breaking down on order the size 
of the smallest grid cell. As a result, errors in the orbits of 
particles may result. However, we believe that our runs are 
short enough to avoid buildup of significant errors. 

In the FLASH runs, to obtain a useful number of star 
particles with long accretion histories, we have taken the 
initial turbulent box and have only run our refinement al¬ 
gorithm (and hence, star particle algorithm) on only one 
octant at a time. This forces us to run eight high resolution 
simulations, each on a difference octant and so allows us 
to treat each octant as a separate distinct simulation. This 
is necessary as FLASH does not have individual timesteps, 
which results in the code grinding to a halt once a single 
region collapses. 


3 RESULTS 

In Figure 1 we show a projection along the z-axis of the en¬ 
tire simulation volume for one of the high resolution octant 
simulations, 2.8 Myr after gravity has been turned on. The 
image shows up to 8 levels of refinement, giving an effective 
resolution of 32768®, or a minimum cell size of 5 x 10“"^ pc. 
Regions that are highly refined are the densest regions, for 
which the image is smoother than the low-density more pixe- 
lated regions. Note that the highly refined regions are limited 
to the lower right, which is the octant that this particular 
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Figure 1. Projection of the density along the z-axis of the entire 
simulation volume. The root grid is 128^ with up to 8 levels of 
refinement, giving an effective resolution of 32 K®. This snapshot 
is taken 2.8 Myr after star formation was turned on. 

simulation focused on. The other seven simulations refine 
the other octants. 

The high density regions are organized into filaments. 
These filaments span most of the simulation box, with 
lengths up to several parsecs and widths of order a few 
tenths of a parsec. Some filaments appear to flow into large 
clumps. This is in accord with many previous simulations, 
e.g., (Padoan et al. 1998; Lee et al. 2015). These clumpy re¬ 
gions have the highest densities and, hence, are the first to 
fulfill the criterion for star particle formation. 

In this section we focus on the regions around two in¬ 
dividual star particles, which we refer to as particle A and 
particle B . 

Particle A formed about a quarter of a parsec away 
from its nearest neighbor star particle. At the end of the run 
it was ~ 736,000 years old and had a mass of ~ 17.5M©, 
although it was still accreting rapidly. 

Particle B formed and remained in isolation. At the end 
of the run, the particle was ~ 512, 000 years old and had a 
mass of ~ 1O.7M0. Throughout the simulation particle B 
had a steady supply of gas. 

3.1 The Run of Infall (vr), Circular (v^), and 

Random Motion (vt) velocities with Radius 
(r) Before Star Particle Formation. 

Figure 2 shows the infall velocity, Ur, circular, v^, and ran¬ 
dom motion, vt, velocities as a function of radius (top panel) 
and the density in a slice of the local volume (bottom panel) 
around the density peak that will form particle A 100,000 
years in the future. In Appendix A, we describe how we 
calculate each of these velocities. 

We will compare vt to what MC15 referred to as a tur¬ 
bulent velocity. Our current definition of vt is simply that 
of a random velocity. We are agnostic about whether or not 


Vt characterizes an isotropic turbulent pressure; close ex¬ 
amination of the velocity field indicates that the random 
motions are not isotropic on the scale of their distance from 
the density peak. It is also clear, however, that vt charac¬ 
terizes a Reynolds stress that does provide a net outward 
support against gravitational collapse. This follows from a 
simple energy argument; the infall velocity in the vicinity of 
the density peak is well below the local free-fall velocity, and 
remains so throughout the simulation, even after a star par¬ 
ticle forms. Thus, some of the potential energy released by 
the infall goes into some channel other than inward motion. 
A fraction of the potential energy release goes into shocks, 
and in our code is effectively removed immediately. At this 
early stage, the rotational motion represents a small fraction 
(< 10%) of the energy at all but the smallest radii. But the 
inward flattening of the green line in Figure 2, and the in¬ 
ward increase seen in later figures, shows that a substantial 
fraction of the potential energy released by the inflow goes 
into random motions. By energy conservation, this fraction 
is not available to the inflow, so that |rir| is smaller than it 
would be if the random motions were not absorbing some 
of the energy. This shows that there is an effective outward 
force on the infalling gas. 

The infall velocity, Ur, and random motion (vt) veloc¬ 
ity are similar in magnitude, and somewhat smaller than 
the Keplerian velocity, vk = \/GM{< r)/r. Note that Ur is 
roughly equal to the sound speed while vt is supersonic. The 
fact that the infall velocity is ~ 25% of the free-fall velocity 
over all radii less than a parsec shows that this system is 
not in hydrostatic equilibrium. The density distribution is 
smooth and filamentary. The run of density versus radius, 
not shown, is a simple power law with a small inner core. 

Figure 3 shows the region around the same density max¬ 
imum some 70,000 years later, 30,000 years before star par¬ 
ticle A forms. Once again the infall velocity is a substantial 
fraction of the Keplerian velocity, showing that the core re¬ 
mains far from hydrostatic equilibrium. However, in the 
innermost regions (inside 0.01 pc) is comparable to both vk 
and Vt, showing that the innermost region is partially rota- 
tionally supported. The density slice, shown in the bottom 
panel of Figure 3, confirms this interpretation, showing a 
disk-like structure with a radius of order ~ 0.01 pc. The 
mass inside this radius is ~ 0.7 Mq. We note that the par¬ 
ticle forms near the tip of a filament (not shown). 

3.2 The Stellar Sphere of Influence 

We begin by developing an operational definition of r*. We 
choose to define r* as the radius where the enclosed mass, 
M{< r,t) is three times the mass of the star, i.e., 

3M4t) ^ M{<r4t),t) (9) 

similar to Murray & Chang (2015). We use the factor of 3 
to ensure that the gravity of the gas dominates the gravity 
of the star.^ In particular, the factor of 3 essentially means 
that the mass in gas is twice the mass of the central mass 

^ The gas in the disk around the protostar is rotationally sup¬ 
ported, so it essentially acts as a part of the star. We include the 
mass of the disk when calculating r* and discuss how we define 
the disk mass in 3.5. 
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Figure 2. The top plot shows the run of velocity with radius 
measured from the density peak; this density peak will develop 
into particle A in 100,000 years. The sound speed is the black 
horizontal line while the infall velocity \ur\ is given by the blue 
triangles, connected by a solid blue line. The green circles con¬ 
nected by a solid green line show vj^ while the black crosses show 
the rotational velocity Vfj^. The red dashed line is the Keplerian 
velocity = y/GM{r)/r. Even at this early stage the structure 
is far from hydrostatic equilibrium, as the infall velocity is ~ 25% 
of the free-fall velocity. The refinement level is Z = 6, which corre¬ 
sponds to a cell size of ~ 2 x 10“^ pc. The bottom plot shows the 
density in a slice along the direction of the angular momentum 
vector centered on that peak. 


(star and disk) and implies that the gravitational accelera¬ 
tion of the gas is twice that of the central mass, which is 
where the dynamical effects of the gas begins to dominate 
the dynamical effects of the central mass. 

Equations (3), (4), and (9) predict that the character of 
the solution should change at r* and that r* increases with 
time. Our numerical results support this prediction. Figure 
4 shows that vt decreases with decreasing radius down to r* 
and then increases with decreasing radius inside the sphere 
of influence. We see that vt reaches a minimum near r = 
The inward decrease in vt{t) is not monotonic near 0.2pc, 
probably due to a shock, as suggested by jumps in both 
the infall and random velocities, and in the density, at r ~ 
0.02 pc. This trend of increasing vt with decreasing radius 
inside r* is repeated in Figure 5. 

We don’t see an increase in the infall velocity for r > r* 
for this object because the star particles are forming about 




Image x (pc) 


Figure 3. The top panel shows the run of velocity around the 
same density peak as that shown in Figure 2, but now only ~ 
30,000 years before the formation of particle A . The color and 
linestyles are the same as in the top panel of Figure 2. The bottom 
panel again shows the density in a slice centered on the density 
peak. The plotted arrows show the velocity in the plane of the 
slice. The longest arrows correspond to roughly 2km s~^. In the 
intervening ~ 70, 000 years since the time shown in Figure 2, an 
accretion disk-like structure has formed, which has a mass of ~ 
0.7 Mq. The radius of the sphere of influence (of the disk) is ~0.02 
pc. All three velocities, |ur|, vx, and increase inward of r*; 
the inflow is disrupted at r ~ 0.015 pc a feature that we interpret 
as a shock, where the flow meets the nascent accretion disk, at 
which point vt also drops in magnitude. At yet smaller radii the 
infall resumes, because at this early time the disk is not yet fully 
rotationally supported. The resolution at the location of the star 
particle has reached the refinement limit Ax = 5 x 10“*^ pc; the 
errors in the calculation of the velocities that are associated with 
the finite resolution are substantial inside r Pi 0.002 pc, so features 
inside this radius are not reliable, and thus not plotted. 


1 pc from the end of a filament, but we do see an increase in 
\ur\ in other particles, see below. 

Comparing Figure 4, which shows the velocity and den¬ 
sity of the same region 24,000 years after star particle A 
forms, with Figure 3 demonstrates that the radius of the 
change in character of the flow associated with r* increases 
over time. In particular, the global minimum of the random 
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Figure 4. The run of velocity (top panel) of particle A 24,000 
years after star particle formation. The color and linestyles are 
the same as in the top panel of Figure 2. This panel shows that 
the disk around the particle is rotationally supported for r < 
5 X 10“^ pc; inside that radius, the black pluses are higher than 
either the green or blue points, i.e., the rotational velocity is larger 
than either the random motionor infall velocity. The radius of 
the sphere of influence is r* ~ 0.06 pc. The star mass is 2.47M0 
and the disk mass is O.SSMq, thus the disk is about 23% the 
mass of the star. The bottom panel is a density slice centered on 
particle A, in the center of mass frame, face-on to the rotationally 
supported disk. The arrows show the velocity in the plane of the 
slice. The longest arrows correspond to nearly 3km s~^. With a 
cell size of Ax = 5 x 10“^^ pc, we have ~ 20 cells across the 
diameter of the disk. 
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Figure 5. The top panel shows the run of velocity for particle 
B ~ 100, 000 years after formation. The radius of the sphere of 
influence is 0.18 pc. The stellar mass is 4.5M0 and the disk mass 
is 2.47M0 so the disk is ~ 55% of the stellar mass. The bottom 
panel shows the density in a slice centered on particle B . 


fraction, say 1/4, of the outer scale. In our simulations, the 
outer scale is given by k = 2, or I//2, and we use solenoidal 
stirring, so that the cascade starts out with no compressive 
component, although one develops as the cascade proceeds. 
In fact we will show in §3.7 that the typical radius of a con¬ 
verging region is more like r « 1 pc in our simulations. 


motion velocity is now at 0.06 pc rather than somewhere 
between 0.01 — 0.02 pc. 

The drop in Ur at large radii in Figure 4 reflects the 
vagaries of the large scale Reynolds stress pressure gradient; 
we already mentioned that this particle is forming near the 
end of a filament. 

Figure 5 shows the velocities and the density in a slice 
centered on particle B, 100,000 years after that star particle 
forms. This star is more isolated than particle A, and as a 
result |wr| increases from r > r* out beyond r ~ 3 pc. This 
is in accord with equations (4) and (5), but it contrasts with 
the result in Figure 4. 

The behavior of \ur{r)\ at large radii is not set by 
the collapse dynamics, but rather by the properties of the 
random motions, most importantly the outer scale of the 
Reynolds stress gradientin particular, we do not expect 
|wr(r)| to be significant on scales larger than some moderate 


3.3 A Fixed Point Attractor for p(r, t) Inside r. 

One of the most striking findings of MC15 was that the run 
of density is independent of time for r < r,. Our simulations 
confirm that finding, as illustrated in Figure 6. The plot 
shows the run of density for two separate times. The dotted 
blue line shows the run of density ~ 40,000 years before 
particle A forms, while the solid green line is the run of 
density ~ 540,000 years after the star particle forms. The 
elapsed time corresponds to nearly two tenths of the mean 
free-fall time of the box, and to many free-fall times at radii 
less than a tenth of a parsec. We emphasize that the density 
can change on the local free-fall time, which is much smaller 
than the global free-fall time (by a factor of 10 or more for 
r < 0.1 pc). We will show that in fact the density inside rd 
does change rather rapidly, after the star particle forms, but 
that for Td < r < rt. the density does not change; see §3.7 
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Figure 6. The run of density for particle A . The dotted blue 
line is the density ~ 40, 000 years before the star forms. The solid 
green line is the run of density Ri 540, 000 years after formation. 
The gap in time corresponds to the free-fall time at a radius of 
0.24 pc. For r < 0.24pc the range of time spanned in the plot 
is more than a local free-fall time, yet the density does not vary 
significantly. There are fluctuation in the density, e.g. the spikes 
around r ~ 0.1 pc In Figure 9 we average over a number of objects 
to remove these fluctuations: we also demonstrate that the density 
in the disc does increase. 

The mean power law slope of the density before the 
star forms (the blue dashed line in the figure) is kp ~ 1.9, 
consistent within the star-to-star variations we see with the 
range kp « 1.6 — 1.8 from equation (3) for r > (since in 
this case r* =0). 

3.4 Mass accretion rate 

In Figure 7 we show the mass accretion rate M as a function 
of r around a star particle (t — t* > 0) and from the corre¬ 
sponding density peak in which the star particle eventually 
formed (t —t* < 0). This plot is taken from a RAMSES sim¬ 
ulation. Before the star particle forms, M decreases inward 
at all radii. 

Following the establishment of the power law solution 
for the density, at t = f*, a star particle forms and the M 
profiles flatten at small radii. An examination of the density 
profile (not plotted) reveals that p oc while for f —f* = 

24kyrs, the gravitational force (and hence Ur) is dominated 
by the central mass for r < 0.1 pc, so that v oc out to 

that radius. We also note that while the M profile is flat, 
it does increase in time as shown by the difference between 
the f — t* = 24.6 kyrs and 154 kyrs curves. All this behavior 
agrees well with the prediction of Equation (5). 

At all times, the accretion rate is either nearly flat or 
increasing with radius, which is a natural result of the near 
balance between gravity and Reynold’s stress support, as 
posited in the theory of MC15. We contrast this with an 
inside-out collapse model, which we exemplify using a Shu 
(1977) solution (blue dashed line) obtained by directly inte¬ 
grating equations (11) and (12) of Shu (1977) at a fixed time. 
The asymptotic behavior of M follows from Shu’s equations 
(15) and (17); recall that x = r/{cst) is a function of the 


Figure 7. The run of M for a star particle in a simulation with 
Nj = 32 cells per jeans length resolution, at a maximum effective 
resolution of 16384^. The lowest solid (green) line is the run of 
M for the density peak that will form the star particle 252,000 
years after the time plotted. The accretion rate is about an order 
of magnitude lower at small radii (say 10“^ pc) than at Ipc. The 
purple line connecting the dots is the run of accretion soon after 
the star particle forms, when the stellar mass is about a solar 
mass. The accretion rate at 1 pc still exceeds that at all smaller 
radii, showing that the collapse is outside-in, not inside-out. As 
an example of an inside-out collapse, we show the accretion rate 
for the Shu (1977) model (the dashed line) for a star of a solar 
mass with A = 5.5. 

radius. In the limit of small x, M approaches a constant. 
However, for large values of x, M{r,t) = —A{A — 2 )cCGx 
(Equation [15] of Shu 1977), i.e., the mass accretion rate falls 
like 1/r at a fixed time at large r as seen in Figure 7. 

In other words, for inside-out collapse models, the accre¬ 
tion rate is monotonically decreasing with increasing radius. 
This is qualitatively different from the prediction of MC15 or 
the results of this work. We note that while we have chosen 
to plot the Shu solution, other collapse solutions (McLaugh¬ 
lin & Pudritz 1997; McKee & Tan 2002, 2003) have the same 
general prohle: the mass accretion rate is roughly indepen¬ 
dent of r at small radii, and decreases with increasing r at 
large radii. 

In summary, at no time do we see any indication of 
an inside-out collapse in our simulated massive star forming 
regions. 

3.5 Rotationally Supported Disks 

Many of the qualitative and even quantitative features pre¬ 
dicted by MC15 are found in our simulations as discussed 
above, including the approach of the density profile inside r* 
to an attractor solution, the minimum in the velocity pro¬ 
file around the sphere of influence, and the expansion of the 
sphere of influence with time. However, our simulations dis¬ 
play additional dynamics that were not modeled by MC15. 

A particularly interesting bit of dynamics neglected by 
MC15 is the development of a rotationally supported disk, 
which we alluded to above. This development is evident in 
the velocity plots, starting from the absence of a disk in Fig- 
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ure 2 to a proto-disk with no central star particle in Figure 
3, to a fairly well developed rotationally dominated disk, at 
r < 0.05 pc in Figure 4. 

We define the outer edge of the accretion disk rd as 
the largest radius where exceeds both \ur\ and vt, that 
is, where the disk is rotationally dominated. The develop¬ 
ment of the disk is best followed by examining the rota¬ 
tional velocity seen in Figures 2, 3 and 4. In the last figure, 
rd ~ 7 X 10“^ pc. We have also used a second definition for 
the disk radius, i.e., where the derivative of the density has 
a sharp drop, see footnote 1. The two definitions of the disk 
radii agree well with each other. 

We note that the disks in our simulation have Vd ~ 
1,000 AU. This is somewhat larger than the radii of the 
largest observed disks, e.g., Padgett et al. (1999) find 
500 AU < rd < 1000 AU. Of course we are simulating mas¬ 
sive star formation, and most observations of disks are of 
nearby, low mass stars. Another factor to keep in mind is 
that we are doing hydrodynamic simulations, so there are no 
magnetic fields, which are believed to be effective at trans¬ 
porting angular momentum; the inclusion of magnetic fields 
might therefore tend to reduce the sizes of the accretion disks 
in our simulations. 


3.6 Gravitationally Unstable Disks 


The plot of M in Figure 7 shows that the accretion rate 
varies little across the transition from the rotationally sup¬ 
ported disk to the radial infall dominated part of the flow 
at slightly larger radii. In other words, the disk is transport¬ 
ing angular momentum efficiently enough so that the disk 
accretion rate matches the rate at larger radii. Since our 
simulations do not include magnetic fields, this efficient disk 
accretion is not due to the magneto-rotational instability 
(Balbus & Hawley 1991, 1998). 

Following Kratter et al. (2010), we suggest that angular 
momentum is transported via gravitational torques. We have 
not yet tried to calculate these torques, but as a first check, 
we have calculated the Toomre Q parameter, as shown in 
Figure 8; recall that 


Q 


TrGrE 


( 10 ) 


In this expression E is the gas surface density of the disk. 
The figure shows Q for the disk around particle B at the time 
shown in Figure 5. For the region 3 x 10“® ^ ^ 6 x 10“^ pc, 

Q is below one, which supports the notion that the efficient 
accretion is due to gravitational torques resulting from a 
gravitational instability in the disk. However, in the next 
section, we find results suggesting that the accretion disks in 
our simulation are not gravitationally unstable at all times. 


3.7 Average Profiles 

Thus far, we have focused our attention on two of our stars 
and shown that their vt, Ur and p profiles are qualita¬ 
tively similar to the profiles predicted in the analytic work 
of MC15. Now we will show that this behavior is generic, 
in the sense that this is true on average over all the star 
particles in our simulations. 

At the end of our base FLASH simulations, we have 



Figure 8. The Toomre Q parameter for particle B , at the same 
time as shown in Figure 5, ~ 100, 000 years after star particle 
formation. From Figure (5) we see that the disk is rotationally 
dominated for r < 2 X 10“^ pc. For 3 X 10“^ $ S 6 X 10“® pc, 
Q < 1. This indicates that the disk is gravitationally unstable at 
these radii, while it is marginally stable at larger radii. Figure 14 
provides a more representative view of disk stability. 


found roughly 60 star particles. To study these systems in 
a generic way, we look at the average velocity and density 
profiles. Motivated by the results of MC15, we average the 
profiles at fixed stellar mass; by fixing the stellar mass, we 
fix r*, and hence the velocity, p, and M profiles. 

For epochs before a star particle forms, it is less clear 
how these profiles should be averaged. However, equation (3) 
predicts that p(r, t) approaches a time independent function 
as soon as any non-pressure supported structure, such as a 
disk, forms. As a result, we elect to follow the methodology 
in Lee et al. (2015) and average profiles at fixed times (10 
and 100 kyrs) before the formation of a star particle. The 
choice of these two times allows us to study the conditions in 
the collapsing region immediately before and well before the 
formation of the star particle, while retaining several (six to 
seven) density peaks and hence reasonable statistics. 

In Figure 9, we plot n as a function of r, 10,000 and 
100,000 years before star particle formation (left plot), and 
for stellar mass M* = 1 and 4 Mq (right plot). The plots 
show that p(r,t) —>■ p{r) for Vd < r < r*, i.e., the density 
approaches the attractor solution, early in the collapse, and 
this profile persists through formation of the star particle 
and well after. This generalizes what we found for our two 
example star particles in §3.3. 

It is important to note that the lack of change in the 
run of density is not due to the fact that we integrate for 
only a few tenths of a global free fall time. To emphasize this 
point, observe that the density a.t r < Vd does increase, while 
the density for Vd < r < does not increase with time. 

The reason for the increase of p(r, t) with time for r < Vd 
is easy to understand: the gas is in a rotationally supported 
disk, which (as Figure 8 shows) is marginally unstable. As 
the central stellar mass grows, the mass of the disk surround¬ 
ing it will grow as well, in such a way that Q ~ Md/M,, 
where Md is the disk mass, is roughly constant. 

The fact that the density inside Vd increases illustrates 


MNRAS 000, 1-22 (2016) 






10 


Murray, Chang, Murray, & Pittman 




Figure 9. Number density as a function of radius at 10,000 and 100,000 years before the star particle forms (left plot) and when the 
star reaches 1 and 4 Mq (right plot). For the left plot, we average over 6 and 7 regions that are within 25% of 10,000 and 100,000 years 
prior to star particle formation. The line corresponding to t — t* = —100,000yrs terminates at r = 0.02 pc because that corresponds 
to the level of refinement for the local density peak (n Ri 10® cm“®) at that time. For the right plot, we average over 14 and 6 regions 
that contain 1 or 4 Mq star particles (within 0.5 Mq). A power law fit to either curve between 0.02 pc and 1 pc gives a power law of 
n oc r~'^p with Kp Ki 1.8. At r = 0.1 pc the free fall time is ^ 250, 000 years, roughly the span of time show across the two panels of the 
plot. For radii between ~ 0.02 pc and 1 pc the density profile does not change appreciably between the left and right plots. The lack 
of change for 0.1 pc < r < 1 pc is unsurprising, since the elapsed time is less than the local dynamical time at these radii. However, the 
same cannot be said about the lack of evolution of p{r,t) in the range < r < 0.1 pc. The collapse theories of (Shu 1977; McLaughlin 
& Pudritz 1997; McKee & Tan 2003) predict that the density should vary with time, but this is not what we see. Note that the density 
profile does increase for r < r^j; compare the green line (the 4Mq profile) versus the blue line (the IMq profile) in the right hand panel. 
This is consistent with mass accreting onto the disk faster than it is transported in towards the star, in such a way that Q Ri 1 at all 
times. The change in the density for r < demonstrates that the density can evolve on the time scale of our simulations. Thus the fact 
that the density does not change for < r < r* between the two plots is striking. It is, however, what is predicted by Equation (3). 


the general point that the relevant time scale for the run of 
density to change can be much shorter than the global dy¬ 
namical time scale. If one had to wait for a global dynamical 
time, the density in the disk would not change over the entire 
course of our simulation, but Figure 9 shows that the density 
in the disk does change over a tenth of the global dynam¬ 
ical (or free-fall) time. Thus the result that p(r, t) —>■ p(r) 
for Td < r < r, is not a result of our short (relative to the 
global dynamical time) integration. 

The one dimensional numerical models in MC15 also 
showed that p{r,t) —^ p(r); MC15 find that the fixed point 
solution is approached from outside-in (see their Figure 1). 
We see the same behavior in the simulations we have run 
with Nj — 16 and with Nj = 32. In those runs we see 
the flattening of the density at small radii and early times, 
before the star particle forms. 

In Figure 10, we show the density probability distribu¬ 
tion function (PDF) of one of our simulations. The black 
line shows the result for the full box. The blue thin dot-dash 
line shows the result when we excise a 1 pc sphere around 
each star particle. Finally, the thin blue dashed line shows 
the PDF of all the 1 pc spheres around each star particle. 
At high densities, the PDF exhibits power law behavior, as 
found by previous workers (Klessen 2000; Kritsuk et al. 2011; 
Lee et al. 2015). Moreover these high density regions are lo¬ 
calized around star particles, as the PDF with 1 pc spheres 
excised around star particles shows (blue dot-dashed line). 
We also note that the regions around star particles are not 
devoid of low density regions, as the PDF of the 1 pc spheres 
around star particles (blue dashed line) shows. Kritsuk et al. 


(2011) first argued that the power law tail of the density 
PDF at high densities is related to the scaling of the density 
with radius; for p oc r““, the density PDF oc For the 

values of a (a « 1.5 — 2) that we expect from analytic the¬ 
ory (MC15) and from previous numerical calculations (Lee 
et al. 2015), we expect the density pdf to scale like p~^ to 
We fit a power law between n — 10"^ — 10® cm~® (red 
dotted line) and find a scaling like in line with these 

expectations. 

The power law shows a break to a flatter slope at 
n « 10® cm“®. A similar break was seen by Kritsuk et al. 
(2011), who argued that at very high densities, the den¬ 
sity PDF flattens due to the presence of disks, which they 
also found. In our simulations, we have found that material 
with n > 10® cm“® always resides within 0.01 pc of a star 
particle. Since 0.01 pc is the typical outer radius of our sim¬ 
ulated disks, this suggests that the highest density material 
is strongly associated with the disk. 

Figure 11 shows the averaged Vr, vt, v^j,, vk as a func¬ 
tion of r before the star particle forms (left) and after (right 
panel). As in Figure 9, we have selected the same fixed times 
(10 and 100 kyrs) before the star particle forms and the same 
fixed masses (1 and 4 M 0 ) after star particle formation. Here 
the dynamics of vt follow quantitatively the behavior of vt 
found in MC15. In MC15, vt scales with radius as at 
very large r, where self-gravity is not important. For exam¬ 
ple, at t = 100 kyrs before the star particle forms, we find 
that vtO ~ r® '^®, in line with Larson’s law, i.e., oc 

However, at t = 10 kyrs before star particle formation, 
one can see that the vt scaling has reversed itself at small 
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Figure 10. The probability distribution function of n. The black 
thin solid line shows the full PDF, the blue thin dot-dashed line 
shows the PDF with 1 pc spheres cut out around star particles, 
and the blue thin dash line shows the PDF within those spheres. 
The red dotted line show the power law oc 

radii (r < 0.1 pc) due to the accumulation of mass in a proto¬ 
disk; the gas in the disk deepens the potential well, but does 
not provide radial pressure support. The figure also shows 
that \ur\ increases inward from r ~ 0.1 pc. 

The reversal of the power-law form of both |t6r| and vt 
as a function of radius tracks the position of r*, as can be 
seen comparing the lines for t = —10 kyrs in the left plot 
with M = 1 and M = 4 M0 in the right plot. This confirms 
another aspect of the MCI5 solution - that as r* moves 
outward with time the inflection point in \ur\ and vt moves 
outward as well, as we found earlier in §3.2. 

The steady outward march of the sphere of influence is 
demonstrated in Figure 12, which shows the run of enclosed 
mass at four different times. At t = —100,000 years, in Fig¬ 
ure 9 the density cusp is not yet in place; correspondingly, 
at small radii the enclosed mass is convex, curving down as 
r decreases. By t = —10,000 years cusp formation is com¬ 
plete, and a small disk has formed, evidenced by a slight 
upward concavity in the enclosed mass profile inside 0.02 pc. 
The radial extent of this upward concavity is increased to 
r « 0.1 pc for M* = 1, and further to r « 0.3 pc by the 
time the stellar mass reaches M* = dM©. The position of 
r* can be inferred by the position in the curves where the 
concave portion of the curve meets the linear portion. The 
concave regions are dominated by a central mass and hence 
are inside of r*. 

The growth of the central mass forces r*(t) outward 
because p{r) is independent of time, and hence the gas mass 
at small radii remains fixed, while Mf{t) grows. 

Returning to the velocities, the fact that vtO ~ 
for t = —100 kyrs in the left plot of Figure 11 shows that 
the turbulence in the initial collapse obeys the same scaling 
law found in non-collapsing regions in the molecular cloud 
(Lee et al. 2015). This suggests that the turbulence in incip¬ 
ient collapsing regions is governed by the same large scale 
turbulent cascade as in non-collapse regions. 

However, the flattening and reversal of vrir) at small 
radii and late times shows that some mechanism other than 


a turbulent cascade is at work at these radii and times. We 
interpret the behavior of vtO as the combined result of 
compressional heating and turbulent decay, as suggested by 
MC15 and by Robertson & Goldreich (2012). 

The relatively large infall velocity demonstrates that, 
even 100,000 years before the proto-disk or star particle 
forms, these regions are not in hydrostatic equilibrium, in 
which Reynold’s stress or turbulent pressure balances the 
force of gravity. This calls into question the assumption 
made by previous analytic models of massive star forma¬ 
tion, such as the turbulent core model. At early times, |ur| 
is between vk/?> and vk, except at r > 1 pc, where the 
clump fades into the ambient molecular cloud. These high 
ratios of \ur\/vK show that hydrostatic equilibrium is not a 
valid description of the star forming regions at any time. 

In fact these plots show that |wr| is of order ux/3 or 
larger at all times for < r < 1 pc. 

At small radii, the fact that p(r, t) = p(r) ~ foj- 

r < r*, combined with the fact that \ur{r,t)\ ~ , en¬ 

sures that M{r,t) = M(t), i.e., the mass accretion rate is 
independent of radius for r < r*. 

This result for the accretion rate was shown previously 
in Figure 7. At early times, {t—t,) ~ —100 Kyrs (the red dot¬ 
ted curve), M decreases by a factor of 20 between r = 0.5 pc 
and r = 0.01 pc because the density profile is still evolv¬ 
ing toward the attractor solution. But for later times M (r) 
is flat at small radii. This demonstrates that the attractor 
solution, once established, imposes a major effect on the ac¬ 
cretion profile. 

While the gas is never hydrostatic, the gradient of the 
Reynolds stress does roughly balance gravity as can be seen 
in Figure 13. The figure shows the rotational support, which 
we define as v^/r (solid blue line), Reynolds stress plus 
thermal pressure support p~^dP/dr = p~^dp{vf /2 -|- (f^ldr 
(dashed line), and total pressure support p~^dP/dr + v’^fr 
(thick red line). We have scaled these quantities to the local 
gravitational acceleration, g = GM{< r)lF. 

Inside of r « 0.01 pc, the gas settles into a rotationally 
supported disk and the support from other sources drops. 
However, the sum of the Reynold stress and rotational sup¬ 
port (thick red line) nearly balances the local gravity. 

The rotationally supported disks in our simulations are, 
on average, roughly marginally stable, as seen in Figure 14, 
which should be compared with Fig. 3 in Kratter et al. 
(2010). For 0.007pc < r < 0.02pc, we find 1 < Q < 1.6. 
Examining individual disks, some of the time the disk is un¬ 
stable and rapidly dumps material toward the central star, 
while at other times the disk is stable, building up mate¬ 
rial to approach marginal stability. For r > 0.02 pc, Q rises, 
though the interpretation of Q as a measure of stability is 
questionable, as the gas is no longer rotationally supported, 
nor is it in a flattened or disk-like configuration. 

3.8 Mass Accretion Rates 

Finally, we discuss the mass accretion rates in our simu¬ 
lation. Previously, Lee et al. (2015) (see also Myers et al. 
2014) found that the star formation efficiency is nonlinear 
in time, with M* oc T. This nonlinear rate is evident in the 
work of previous workers, but was often interpreted as an 
initial transient (Padoan & Nordlund 2011). MC15 showed 
that M*(t) ~ is a natural consequence of the density ap- 
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Figure 11. v^, vj-, and as a function of r at 10,000 (thin lines) and 100,000 (thick lines) years before the star particle forms (left 
plot) and when the star reaches 1 and 4 Mq (right plot). The averages are over the same regions as those used in producing Figure 9. 
The infall |ur| and random vx velocities show the behavior predicted by the theory of adiabatic turbulent heating for times later than 
— 10,000 years: at large radii, where |Mr| is small, vx > |iir| and vx decreases inward, but more slowly than in non-collapsing regions; 
p Ri 0.2 rather than p = 0.5. Inside r,, where \ur\ > vx (or |lf| > vx in the notation of Robertson & Goldreich (2012)) vx increases 
towards |ur| as r decreases, with both increasing inward. 






Figure 12. Mass of gas and stars as a function of r at 10,000 (thin lines) and 100,000 (thick lines) years before the star particle forms 
(left plot), and when the star reaches 1 and 4 Mq (right plot). Averages are as described in Figure 9. 


proaching an attractor solution and the scaling of the infall 
velocity with the Keplerian velocity at small radius, as we 
have clearly demonstrated in this work. 

First we address the question of whether the M* oc 
phase is an initial “transient”. Tackenberg et al. (2012) and 
Traficante et al. (2015) estimate the lifetimes of massive star¬ 
forming clumps found using the Apex telescope and the Her- 
schel telescope, respectively. Tackenberg et al. (2012) iden¬ 
tify clumps with column densities Eg > O.lgcm”^, masses 
up to 10® Mq. Since they have a fairly complete catalog 
of such clumps, they can estimate the typical lifetime of a 
clump by comparing to the number of massive stars formed 
in the Milky Way every year. They find a mean clump life¬ 
time of 6 X 10"^ yr, and a clump free-fall time of « 1.5 x 10® yr; 
the clumps live 1 free-fall time. 

Similarly, Traficante et al. (2015) identify clumps with 
sizes ranging from 0.1 — Ipc, masses ranging up to IO'^Mq. 


They estimate an upper limit lifetime for the starless phase 
of 10® yr for clumps with M > 500Mq, and a ratio of starless 
to total clumps (the rest of the clumps host protostars) of 
39%. Thus the total lifetime of clumps in their mass range 
(above 500Mq) is ~ 2 — 4 x 10® yr. The clumps in their 
sample have 10'*cm“® uh < 10 ®cm“®, so free-fall times 
1.6 X 10® yr < th < 5 x 10® yr. The clumps live 0.2 — 2.5 
free-fall times, similar to the estimate of Tackenberg et al. 
(2012). 

Thus the lifetimes of massive star forming clumps, when 
measured in units of free-fall times, is similar to the lifetimes 
of GMCs, again measured in free-fall times, e.g.. Blitz et al. 
(2007), who find that GMCs live 2-3 free-fall times. 

Our simulations run for only a fraction of a free-fall 
time, but those of other workers have often run for two to 
three (Wang et al. 2010; Padoan & Nordlund 2011), or, in 
some cases, up to five free-fall times. The simulations are of- 
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Figure 13. The ratio of rotational (solid line) support, Reynolds 
stress and gas (dashed line) pressure support to the local gravi¬ 
tational acceleration g = GM{< r)lr’^, as a function of r when 
the star reaches 4 Mq. 



Figure 14. The average (over many disks) of Toomre Q as a 
function of r when M* Ri 4 Mq. 

ten halted when ~ 10 — 30% of the mass has been converted 
into stars, since that is a rough observational estimate of the 
maximum fraction of clump gas turned into stars, e.g., Lada 
& Lada (2003). In most cases this star formation efficiency 
is reached in one or two free-fall times, while the M*(t) ~ T 
behavior is still apparent from the plots. In the case of Feder- 
rath (2015), in the MHD run with stellar wind/jet feedback, 
the M{t) ~ behavior ceases after 4 free-fall times, when 
the star formation efficiency is about 15%. Since it is un¬ 
likely that star forming clumps live so long, the simulated 
star formation rate is probably too low. 

We conclude that the time scale over which the T scal¬ 
ing is seen in simulations is similar to the lifetimes of massive 
star forming regions, so that, while the behavior we focus on 
may be of short duration, it is not “transient”. 

In Figure 15 we show the total M, as a function of time 
since the first star particle was formed, t*. This is exactly the 


same analysis as Lee et al. (2015). However, because the sim¬ 
ulations are distributed among the eight different octants, 
each with a different star formation time, we produce the 
total SFE history as follows. First, we analyze the simula¬ 
tions to find the earliest time at which a star particle formed, 
which we define as t,. We then look at all the simulation to 
find the earliest time at which a simulation ended or tend, 
which defines the time over which all our simulations have 
data. Because each snapshot for each simulations are taken 
at different times, we define a number of times at fixed in¬ 
tervals between t* and tend and interpolate the total stellar 
masses for each simulation on those times. These masses are 
then summed to produced M*(t), which we plot in Figure 
15 

As shown in Figure 15, M*(t) grows roughly linearly for 
~ 100, 000 yrs after the first star particle is formed. However 
once the total stellar mass reaches about M* > IOMq, at a 
time t — t* >100 kyrs, M* oc t^. At this stage, M*/Mgmc ~ 
0.001. This agrees well with the results of Lee et al. (2015), 
who found M* ~ T for stellar masses between M*/Mgmc ~ 
0.015 and 0.3. Due to the computationally expensive nature 
of our much higher resolution simulations, even given the use 
of AMR, we are not able push our simulation to the same 
total M,/Mgmc as Lee et al. (2015) were able to in their 
fixed grid, but much lower resolution, simulations. However, 
our simulations do show that their simulations were already 
at sufficiently high spatial resolution to recover the scaling 
relation. 

The reason for this is not hard to find. Figure 9 shows 
that at r « 1 pc the density has already settled onto its 
time-independent form, while Figure 11 shows that the in¬ 
fall velocity is scaling as for r < 0.3 — 0.7 pc, with the 

smaller value corresponding to the time of star particle for¬ 
mation, and the larger value to times for which M* « IMq. 
As long as a simulation resolves this radius (which corre¬ 
sponds roughly to r*), it will recover the M* ~ T scaling. 

The slower growth at earlier times is due to fact that, 
at the time of star formation, the infall velocity is non-zero, 
despite there being no star to attract the gas; see Figure 
11. In other words, the initial infall velocity is larger than 
it takes time before the Reynold stress can slow 
the infall to the steady state value given by equation (4), 
and hence before the mass accretion rate settles onto the 
steady state value given by equation (5). 

The same comments apply to the accretion rates of in¬ 
dividual stars; individual stars start out accreting mass at a 
roughly constant rate. At later times, when they have sub¬ 
stantial masses, the accretion rate grows linearly in time. 
This happens only with the most massive star particles in 
our simulation. 

The total number of star particles also grows roughly 
linearly with time; the combination of this linear growth 
in number of stars, together with the roughly linear mass 
growth of most of the stars, produces the over all M* ~ 
scaling we see for the simulation region as a whole. In regions 
where the summed accretion rate is highest, individual star 
particles are not able to accrete all the collapsing mass, lead¬ 
ing to the formation of new star particles in the immediate 
vicinity. In other words, our simulations produce clustered 
star formation. This is similar to what has been found in 
other recent simulations (Lee et al. 2015; Gong & Ostriker 
2015). 
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Figure 15. The total mass in stars plotted as a function of time 
since the first star particle formed, i*. The final value of the total 
stellar mass M* Ri ISOMq is about 1% of the total gas mass in 
the box, while the final t — t* Ri 0.6 Myrs, about 15% of the free 
fall time 3.8 Myrs for the mean density of the box. The dashed 
green line shows a slope of 2. 

We interpret Equations (3-5) as a description of star 
cluster formation; the total mass of the cluster will grow as 
the most massive stars may spend a significant fraction 
of their accretion history growing as , but many of the less 
massive stars will not undergo snch rapid growth in their 
accretion rate. 

4 DISCUSSION 

4.1 Basic Results of this Work 

We begin with a summary of our results. 

4-1.1 Collapse is not self-similar 

First, in our isothermal, driven turbulence simulations, star 
formation is not a self-similar process. Two length scales, 
in addition to the radius of the outer boundary condition 
or turbulence outer scale, and the stellar radius, enter the 
problem, the Keplerian or outer disk radius r^, i.e., the ra¬ 
dius at which the gas becomes rotationally supported, and 
the radius r* of the stellar sphere of influence, the smallest 
radius at which the gravity of the gas dominates the grav¬ 
ity from the star and disk. The existence and significance of 
rd has been known since the time of Kant and Laplace; the 
recognition that r* plays a role in star formation is recent, 
so we concentrate on the effect of r* on the dynamics in 
what follows. The value of r* increases monotonically with 
time, since the stellar mass increases monotonically, while 
rd may vary in or out with time, depending on the (turbu- 
lently determined) distribution of angular momentum of the 
accreting gas. 

The non-self-similar behavior of the collapse is most 
strongly reflected in the variation of vt with radius; for 
rd < r < r,, the random motion velocity is a decreasing 
function of radius vrir) ~ r^ with p « —1/2, while for 


r* < r, it is an increasing function of r, scaling like r°'^. 
Similarly, the infall velocity |rir| ~ r~^^^ for rd < r < r*, 
while for r* < r it is fiat or even increasing outward (as in 
the top panel of Figure 5), with substantial variations both 
from particle to particle and at different times for the same 
particle, due to the vagaries of the turbulent flow at large 
radii. 


4-1.2 Density approaches an attractor solution 

Second, we find that inside the sphere of influence of the star, 
the density remains constant over several to tens or even 
hundreds of (local) dynamical or infall times; for rd < r < 
r*, p(r, t) —> p{r). This is illustrated by Figures 6 and 9. One 
implication of this result is that one cannot use observations 
of the free-fall or crossing time of collapsing structures to 
infer either the age or lifetime of those structures. 

The fact that p{r,t) = p(r) ~ ^-“3/2 fgj- ^ com¬ 
bined with the fact that \ur{r,t)\ ~ ensures that 

M{r,t) = M{t), i.e., the mass accretion rate is independent 
of radius for r < r^ (see Figure 7). 

Since \ur\ increases with stellar mass and hence with 
time where r < r* (see the second panel in Figure 11), 
while p is fixed, M{t) increases with time, a result seen in 
many previous papers (Padoan & Nordlund 2011; Bate 2012; 
Krumholz et al. 2012b; Federrath & Klessen 2012; Myers 
et al. 2014), although this fact was usually not commented 
on. After some initial transient behavior, we find M{t) ~ t^ 
(Figure 15) in line with the results of Lee et al. (2015). 

4-1.3 Partitioning of the Collapsing Region’s Potential 
Energy 

Our third result is to show how the potential energy released 
in collapse is partitioned. In our simulations, the support 
from random motions slows the rate of infall, so that |Mr| 
is significantly smaller then the free-fall velocity, but large 
enough to maintain vt at a sufficient level that the accelera¬ 
tion due to the Reynolds stress is close to the acceleration of 
gravity (Figure 13). In contrast, Sur et al. (2010) and Fed¬ 
errath et al. (2011) find an infall velocity which is equal to 
the free-fall velocity just outside their core. The dynamics 
in their simulation is very different to the dynamics in ours; 
in their case the acceleration due to the pressure gradient is 
negligible compared to the acceleration due to gravity. Their 
initial conditions incorporate transonic turbulence, but no 
driving. Since the collapse in their simulation does not take 
place for roughly four free-fall times, by the time the collapse 
starts, the turbulence is subsonic. 

Because the infall in our simulations is typically super¬ 
sonic, some of the kinetic energy can then be converted into 
thermal energy by shocks; of course, even in the absence of 
shocks, normal (molecular) adiabatic heating will convert a 
small fraction of the liberated potential energy into heat. 

Figures 3, 4, and 5 show that for r > r* the bulk of 
the potential energy goes into random motion, and thence 
into shocks. Inside r*, but outside of the disk, the potential 
energy that is not immediately radiated is shared roughly 
equally between the infall and random motions. Inside of 
the disk, the potential energy is converted to rotational and 
random motion, in roughly equal measure, and thence to 
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thermal emission. At any radius, the ratio of kinetic energy 
to potential energy is typically around a quarter to a half, al¬ 
though at large radii (r > 1 pc) the turbulent kinetic energy 
can exceed the potential energy: on the scale of our box, 
the ratio vt/C GMgmc /roMC « cdjf « 1 , where Mgmc 
is the total mass contained in the simulation volume and 
roMC = 8 pc is the ’’radius” of the simulation volume, i.e., 
half of the side of the box. This ratio scales as ~ 1/v^i so 
at r « 1 pc the ratio is ~ \/ 8 - 

We also found evidence that the ratio of radial to 
transverse random motion depends on whether the collapse 
leads to radial compression (for r > r.) or dilation (for 
r < r*), and on the tendency for hydrodynamic turbulence 
to isotropize motions down the cascade; see Figure C2 and 
the discussion in §C1. 

4-1-4 Modification of Larson’s Law 

Our fourth result is the confirmation that the adiabatic heat¬ 
ing of the turbulence alters Larson’s law. On large scales or 
away from collapsing regions, Larson’s law is vt ~ r^ with 
p « 0.4 — 0.5. It emerges naturally from the decay of su¬ 
personic turbulence that is driven on large scales. We find, 
as did Lee et al. (2015), that in rapidly collapsing regions 
the decay of vt ~ r^ with decreasing radius is slowed for 
^ ^ 1PC- Least squares fits to VtO for over this range 
of radii result in exponents between 0.1 < p < 0.4, with 
an average around p « 0 . 2 , in fair agreement with the pre¬ 
diction of equation (4). Inside of r*, we find p = —1/2, as 
predicted by equation (4), representing a reversal of Larson’s 
law. 

4-1.5 Collapse does not proceed in an inside-out manner 

A hfth result is that the gathering and accretion of mass 
starts from large scales, and that, both before and after a 
star particle forms, M{r,t) is larger at large r than it is at 
small r; in other words, the collapse proceeds in an outside- 
in manner. The first point, that the accretion starts from 
large scales, is illustrated by Figures 2, 3, and the left panel 
of 11 , which show that |'Ur(r)| ~ (l/3)tiK(r) out to r ~ Ipc 
or farther, and is trans- or supersonic tens or hundreds of 
thousands of years before the density cusp forms, and hence 
before star or even disk formation starts^. 

Figure 7 shows that just after the cusp/star particle 
forms, the mass accretion rate is actually larger at larger 
radii. This behavior is the opposite of that predicted by 
inside-out collapse models, either that of Shu or of the tur¬ 
bulent core model. In Figure 7, this inside-out behavior is 
illustrated with a Shu-type solution, shown by the dashed 
line; in that solution, the mass accretion rate decreases with 
increasing radius, in contrast to the results of our simulation. 

Figures 5 and 11 (right panel) show that just after and 
well after the star forms, the surrounding region is also far 
from hydrostatic equilibrium. 

We conclude that there is no indication of gas in hy¬ 
drostatic equilibrium prior to, during, or after star particle 

^ We remind the reader that the scale of the infall region in our 
simulations, and possibly in the ISM of galaxies, is a simply a 
fraction of the driving scale of the turbulence. 


formation in our simulations. Nor is there any indication of 
inside-out collapse. 

The violation of self-similarity and evidence against 
inside-out collapse shows that the assumptions made by pre¬ 
vious analytic collapse models (Shu 1977; Myers & Fuller 
1992; McLaughlin & Pudritz 1997; McKee & Tan 2003), are 
not fulhlled in our simulations. In addition, the collapsing 
regions in our simulations do not start from a hydrostatic 
equilibrium. 

4.1.6 The magnitude of the pressure gradient term is 
comparable to that of the gravity term for r > rd 

Figure 13 shows that the acceleration due to the pressure 
gradient is comparable to the acceleration of gravity for < 
r < r*. Thus Reynold’s stresses slow the infall compared 
to the free-fall rate, i.e., \ur{r,f)\ < \JGM(r)/r. At small 
radii (r < Vd), the rotational support becomes important 
and the support from vt becomes much smaller than the 
radial component of gravity, but comparable to the vertical 
component in the disk. 

4-1.7 The total stellar mass increases as T 

The total stellar mass in our simulation region, and in in¬ 
dividual star forming sub-regions, increases as the square 
of time after the first star (in the box, or in the individual 
star forming region) forms. Low mass (less than a few so¬ 
lar masses) stars have M*(t) ~ C, with 0 < 7 < 2 , with a 
typical value 7 ~ 1 , but the total number of low mass stars 
N{t) ~ t, so that the total mass in low mass stars grows as 
T. High mass stars, which tend to sit at density peaks (or 
at the bottom of potential wells) have M*(t) ~ t^. 

4.2 Comparisons to Observations 

Caselli & Myers (1995) showed that massive star forming re¬ 
gions have shallower line width-size relations than the clas¬ 
sical Larson result, i.e., vtO ~ r^ with p = 0.21 ± 0.03, 
compared to p « 0.53 ± 0.07 in low mass star forming 
regions. Plume et al. (1997) also found that Larson’s law 
breaks down in massive star forming regions, i.e., their mea¬ 
sured line widths are larger for a given source size than those 
found in low mass star forming regions. As noted in §4.1.4, 
we find the same behavior in our simulations, and we inter¬ 
pret this as the effect of adiabatic heating in a collapsing 
flow at r > r,. 

In addition. Plume et al. (1997) plotted the mean veloc¬ 
ity dispersion as a function of number density, which they 
derived from an excitation analysis of CO. They found that, 
contrary to expectations, the velocity dispersion increased 
with increasing density, which is opposite to the expectation 
based on Larson’s law or supersonic turbulence driven from 
large scales. They concluded that the conditions in dense 
star forming cores are different from the rest of the cloud. 
The simulations presented here, and the analytic results of 
MC15, show the same behavior. In particular, the theory 
suggests that the enhanced turbulence or velocity dispersion 
at small radii in dense star forming regions is the result of 
gravitational collapse adiabatically heating the turbulence. 

We hnd qualitative agreement between the observations 
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of Plume et al. (1997) and our results, i.e., enhanced line- 
widths at high densities, which are associated with smaller 
radii. Performing a more detailed comparison is more diffi¬ 
cult as we have selected regions with the same stellar mass, 
whereas the stellar mass in Plume et al. (1997) is not well 
known. However, it is promising that the linewidths in our 
simulations are of similar magnitude and show the same 
trend with density as do the observations. 

There are now numerous measurements of infall at large 
radii ~ 0.1 — 1 pc in the literature. For example, Csengeri 
et al. (2011) Cygnus X, D = 1.7kpc see infall \ur\ ~ 0.1 — 
0.6kms“^, ot ~ 0.6 — 2kms“^ at r « O.lpc, n ~ 10® — 
10®. Other examples include Ragan et al. (2012, 2015) and 
Peretto et al. (2013). 

Infall is also seen on larger scales, r « 1 pc, by Wyrowski 
et al. (2016), who observe |ur.| in the range 0.3 — 3kms“^, 
corresponding to a fraction of the free-fall velocity (1.4 the 
Keplerian velocity) of 0.03—0.3. In words, the gas at r = 1 pc 
is not in hydrostatic equilibrium, nor is it in free-fall. The 
turbulent velocity in the same clumps at the same radii is 
comparable or slightly in excess of the infall velocity, ut ~ 
1.0 - 2.3kms“®. 

Ho & Haschick (1986), Klaassen & Wilson (2008), and 
Klaassen et al. (2011) see infall at three different radii, r « 
0.5 pc, r « 0.3 pc, and r « 0.03 pc using different molecular 
tracers in the same object, GlO.6-0.4. The infall velocity is 
large at r, small at r « 0.3 pc, and large again at r « 0.03 pc. 
As noted by MC15, this is in qualitative agreement with the 
picture of adiabatically heated turbulence. 

4.3 Missing physics 

Our current understanding of star formation suggests that 
the effects of magnetic helds, radiative and proto-stellar out¬ 
flow feedback from stars, and the equation of state of the gas 
can all have signihcant effects on both the rate of star forma¬ 
tion and the initial mass function (IMF) of the stars. We do 
not include any of this physics in the simulations described 
in this paper. 

It is often argued that the turnover in the IMF, some¬ 
where between 0.2 and 0.6 M©, is associated with the ther¬ 
mal state of the gas in the collapsing region. If so, then our 
use of an isothermal equation of state suggests that the IMF 
found in our simulations is likely to be in error, so we have 
not discussed our computed IMF. However, as Figures 2- 
5 and 11 show, both lu^l and vt exceed Ca, except at the 
earliest times (~ 100, 000 yrs before a star forms), and then 
only for r < O.lpc, so that the gas pressure does not domi¬ 
nate the dynamics in most regions and most of the time. Of 
course we do include the effects of gas pressure, so even in 
those regions and those times, our simulations capture the 
dynamical effects to lowest order, aside from, as we have just 
said, from fragmentation effects on the smallest scales. 

We have undertaken and made some preliminary analy¬ 
ses of magnetohydrodynamic simulations, which we will re¬ 
port on in future publications; as seen by other authors, we 
hud that magnetic helds slow the star formation rate. But 
the runs of density and velocity have the same qualitative 
form in our MHD simulations as in the hydro runs presented 
here, and the MHD runs also give M*(t) ~ 

Like magnetic helds, feedback from protostellar outhows 
are seen to slow the rate of star formation, e.g., Wang et al. 


(2010); Federrath (2015). But those authors also hnd that 
dMt/dt increases with time even in runs that include out¬ 
hows. 

Radiative feedback will also affect both the IMF and, 
for massive enough stars, the dynamics of the collapse at 
late times (after massive stars have formed). 

All the hgures we show present results for stars with 
masses no larger than about AMq. To estimate the effects 
of radiation, we compare the force from the Reynold stress 
Ft = ‘i'KppVT, to the radiation force Lie. From Figure 11, 
the (averaged over many stars) vt is slightly in excess of 
lkms“^ at r = 0.01 pc, while from any of the density hg¬ 
ures the density is p « 5 x 10“^®g/cm®. The force from 
Reynold stress is then f « 4 x 10^® dynes. The luminos¬ 
ity of a 4 solar mass star on the zero age main sequence is 
L « 2 X 10®®ergs“^ (Schaller et al. 1992), so the radiation 
force Ljc « 3 x 10^® dynes, about a 10% effect. The force 
from Reynolds stress increases outward, see Figure 13, so 
this statement holds at larger radii as well. 

Thus we expect that the effects of radiation pressure are 
not particularly signihcant in the situations we report; the 
run of density and infall velocity, and hence the M„{t) ~ 
scaling should not be affected, at least up to the times 
we are reporting on. We note, however, that this estimate 
neglects the effect of radiative or ionization heating which is 
an important feedback mechanism. 

Simulations including radiative feedback support this 
simple analysis. Figure 15 of Myers et al. (2014) shows that 
in their simulations, which include feedback from both pro¬ 
tostellar outhows and radiation (as well as magnetic helds), 
the stellar mass increases as the square of the time, up to 
masses of 4.5 solar masses. Earlier work by the Berke¬ 
ley group found similar results, forming stars with 10 solar 
masses, with M„{t) even for such massive stars, see Figure 
13 of Krumholz et al. (2012b). Their simulations included 
radiative ehects, but no proto-stellar winds. 


5 CONCLUSIONS 

Motivated by recent analytic (MC15) and numerical 
(Padoan et al. 2014; Lee et al. 2015) results, we perform deep 
AMR simulations of star formation in self-gravitating con¬ 
tinuously driven hydrodynamic turbulence. We show that 
two length scales emerge from the process of star forma¬ 
tion, r* and Vd, and demonstrate that these length scales 
are clearly associated with physical effects. In particular, the 
character of the solution changes at ra.{t), inside of which 
(but outside Vd) |ur-| and vt are both oc outside of 

r*, Vt ~ (with p « 0.2), while |ur.| is on average about 
constant. We emphasize that the length scales at which the 
character of the solution changes are time dependent. As the 
star grows in mass, the radius where the stars’ gravity ex¬ 
ceeds the gravity of the surrounding gas increases outwards 
away from the star, r*(t) oc M*^®(t). The disk radius, Vd, 
also changes as a function of time as a result of the advec- 
tion and transport of angular momentum from largescales 
to small scales (and vice versa). 

We also found that the density prohle evolves to a hxed 
attractor, p(r, t) —>■ p(r) in line with the results of MC15 and 
the earlier numerical results of Lee et al. (2015). 

Our results strongly support the basic premise of MC15, 
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that turbulence is a dynamic variable which is driven by adi¬ 
abatic compression (Robertson & Goldreich 2012), and that 
the turbulence in turn acts to slow the collapse. We note, as 
did MC15, that observations of massive star forming regions 
also find vt oc with p ~ 0.2 — 0.3, and that at small radii 
or high density, vt increases with increasing density, as seen 
in observations of massive star forming regions (Plume et al. 
1997). We hnd these departures from Larson’s law only in 
collapsing regions in our simulations. We also show that the 
acceleration due to the pressure gradient is comparable to 
that due to gravity at all r > r^. As a result, the infall ve¬ 
locity is substantially smaller than the free fall velocity even 
very close to the star or accretion disk. Inside r^j, rotational 
support takes over and as a result Ur and vt both decrease. 

Our simulations capture rotational dynamics that 
MC15 did not capture in their 1-D model. In particular, 
we hnd the development of rotationally support disks at 
rd ~ 0.01 pc. These disks have radii comparable to or 
slightly larger than disks seen around young stars in Taurus 
(Padgett et al. 1999) in which stellar feedback effects are 
minimal, and where the undisturbed disks are larger than 
in more active star forming regions such as Orion, where 
the disk radii are ~ 100 AU (Williams & Cieza 2011). This 
is despite the fact that we do not include magnetohydro¬ 
dynamic effects in our numerical computations; large scale 
magnetic helds may transfer angular momentum away from 
these disks, shrinking them. 

Like the disks modeled by Kratter et al. (2010), our sim¬ 
ulated disks are marginally gravitationally stable, suggest¬ 
ing that large scale gravitational torques are responsible for 
transport of material and angular momentum in our simula¬ 
tions; this may also be true at early times in real protostellar 
disks. 

We have shown that the assumptions made by previous 
analytic collapse models (Shu 1977; Myers & Fuller 1992; 
McLaughlin & Pudritz 1997; McKee & Tan 2003), are not 
fulfilled in our simulations. In particular, the collapsing re¬ 
gions in our simulations do not start from a hydrostatic 
equilibrium, nor do they show any evidence of inside-out 
collapse. The gathering of material before collapse, i.e., be¬ 
fore the central cusp in the density power law is formed, 
involves transonic bulk motions and supersonic random mo¬ 
tions (see Figure 11). The accretion of mass starts at large 
scales (r ~ 1 pc) with large initial infall velocities. In addi¬ 
tion, we find that vt scales differently in collapsing regions 
as opposed to the rest of the simulation box, whereas the tur¬ 
bulent collapse models (McLaughlin & Pudritz 1997; McKee 
& Tan 2003) assume that the scaling of vt with r remains 
fixed. 

Finally, we close with a brief discussion of how our re¬ 
sults relate to turbulence regulated theories of star forma¬ 
tion. Here we find several points of disagreement. First, we 
find that the star particles accrete continuously from the sur¬ 
rounding large scale turbulent flow; there is no hydrostatic 
“core” that is cut-off from the turbulent medium. Second, the 
density distribution does not remain log-normal, but rather 
develops a power law tail that is directly related to the den¬ 
sity profile (Kritsuk et al. 2011; Lee et al. 2015). Third, the 
fact that the density profile approaches an attractor solu¬ 
tion that scales like for r < r* and Ur scales with 

the Keplerian velocity guarantees that M is constant with 
radius and M* oc t and hence a non-linear star formation 


efficiency, i.e., M, oc T results. This is in contrast with tur¬ 
bulence regulated theories of star formation that predict a 
constant star formation rate, i.e., M* = const and, hence, a 
linear star formation efficiency M* oc t. 
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APPENDIX A: CALCULATING THE RANDOM 
MOTION AND ROTATIONAL VELOCITY 

In this appendix, we discuss how we calculate the random 
motion, vt, infall, Ur, and rotational velocities, V 4 , from 
the full three dimensional numerical solution. To begin, we 
adopt a series of concentric, logarithmically spaced, spherical 
shells around either a star particle (if available) or around 
a density maximum in the case where the star particle has 
not yet formed. We then removed the bulk velocity from 
these shells by first calculating the enclosed mass, M(< r), 
and momentum in each sphere P(< r), then dividing the 
two to find the bulk velocity, V = P(< r)/M{< r) of the 
sphere of matter. We then subtract this bulk velocity from 
the corresponding shell. 

We also tried defining the bulk velocity using the total 
momentum in each of the spherical shells (rather than in the 
enclosed spheres), and found very similar results. 

We then subtract the bulk velocity from the raw velocity 
of each cell in the spherical shell. We denote the result by v. 

Having removed the bulk velocity, we then calculate the 
radial infall velocity, Vr = v • f per cell, where v is the 
velocity of the gas in a cell and f is the radial unit vector 
(with the origin at the location of the star or local density 
peak). Finally we find Ur =< Vr > as the average over the 
spherical shell, where 

(Ur) = ^i’eeii'Vr^i I Mshell (-^1) 

denotes a mass weighted average over each spherical shell, 
and Maheii is the mass of the shell. The sum is over all the 
cells in the thin spherical shell. 

To calculate the velocity in the (j> direction, where 4> is 
defined by taking the a axis along the angular momentum 
vector of the shell, we first calculate the angular momentum 
Lsheii = /shell ^ ^ 'vdm, where m is the mass in a cell. We next 
calculate the moment of inertia tensor I of each spherical 
shell. In component form, I is 

lij = / {SijP — XiXj)dm (A2) 

J shell 

We then find the rotation vector Cl by inverting 

L = ID, (A3) 

e.g., (McKee & Zweibel 1992). Next we calculate the rota¬ 
tional velocity in each cell from 

= D X r. (A4) 

This amounts to assuming that the gas in each spherical shell 
rotates rigidly; in other words we are averaging over the ran¬ 
dom motions in the shell. Finally, we calculate the spherical 
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shell average as =< |v^| >, i.e., the mass weighted aver¬ 
age of the norm of v,^. 

Armed with the coherent infall (ur) and rotational (u^) 
velocities, we define the remaining velocity as the random 
motion velocity (per cell) as 

Vt = V — UrV — V 0 (A5) 

and the spherical average as vt =< |vt| >. As a check that 
we were accounting for all of the velocities, we added the 
velocities in quadrature: Usum = + v'^ and veri¬ 

fied that it traces the mass weighted average total velocity, 
Vtot =< |v| > accurately. 

APPENDIX B: FILAMENTARY OR 
SPHERICAL ACCRETION? 

Figures 2-5 show that the density in the vicinity of collapsing 
regions is decidedly non-spherical. Despite this, the results of 
MC15 appear to describe the accretion process well. For ex¬ 
ample, in those same Figures we have shown mass-weighted 
infall, random motion, and rotational velocities, the first two 
of which behave as predicted by MCI5 (they made no pre¬ 
dictions for Vtf,). 

To understand this better, we examine how M depends 
on p, and how both are distributed on the sky as seen by 
the accreting particles. The left plot of Figure B1 shows a 
histogram of cumulative M{p)/Mtot through two spherical 
shells at r = 0.5 pc and 0.05 pc, as a function of p/(p), where 
(p) denotes the density average over the (finite thickness) 
shell. We show average histograms when the central star 
has a mass M* = IMq (dashed lines) and M, = 4 M 0 (solid 
lines). The plot shows that 50% of the accretion through the 
sphere occurs via gas that has a density less than 2-5 times 
the average density of the shell, where the low end of this 
range occurs at small radii at late times, with the high end 
occurring at large radii and early times. 

Since the mean density at r = 0.5 pc is (p) « 3 x 
10“^^gcm“^, see Figure 6 or Figure 9, an examination of 
Figure 5, where gas with three times the mean density is 
depicted by dark green (and less dense gas is blue), shows 
that more than half of the accretion is coming from gas that 
covers most of the sky as seen from each of those accret¬ 
ing particles; most of each slice is colored blue. If we take 
filaments to consist of gas that is colored light green or yel¬ 
low (with p > or ~ 3 times the mean density (p) 

inside r = 0.5pc), the filaments account for less than half 
the accretion. 

A similar statement holds for the accretion inside r = 
0.05 pc, shown as the thin lines in Figure Bl. 

To see more quantitatively how this gas is distributed 
on the sky, we plot in the right panel of Figure Bl the cumu¬ 
lative solid angle as a function of p/(p), again for r = 0.5 pc 
(thick lines) and for r = 0.05pc (thin lines). Roughly 90% 
of the sky is covered by gas that is at three times the av¬ 
erage shell density or lower, consistent with the qualitative 
analysis in the previous paragraph. 

Figure B2 shows a histogram of the cumulative nor¬ 
malized M as a function of the cumulative normalized solid 
angle. The plot shows that half the accretion occurs over 
about 10% of the sky where the density is ~ 3 or more times 
the mean density of the spherical shell. So, while about half 


the gas accretes from over most of the sky, and at about the 
mean density, very dense gas entering the sphere from a very 
small covering fraction of the sky contributes the other half 
of the total accretion budget. 

Thus, while the filaments are readily identifiable by eye, 
and are important sources of accreting gas, much of the ac¬ 
cretion (and much of the mass) lies in gas that is more nearly 
spherically distributed. 


APPENDIX C: STAR FORMATION CRITERIA 

The majority of our simulations used a simple density con¬ 
dition of three times the Truelove condition (Equation 8) 
at the maximum refinement level inspired by the sink parti¬ 
cle formation criteria of Padoan & Nordlund (2011) as dis¬ 
cussed in §2. We have experimented with additionally in¬ 
cluding the sink particle checks of Federrath et al. (2010b) 
to check the robustness of our results to these additional 
checks. In Figure Cl we show the run of velocity in a simula¬ 
tion in which we included the star particle formation checks 
used in the default used in FLASH. The results do not dif¬ 
fer significantly from runs lacking such checks. For example, 
both show that the stellar mass increases like t — t* squared, 
M*(t — t*) oc {t — . There are however stochastic varia¬ 

tions in the stellar mass ratio from runs with and without 
the extra checks. The mass ratio at a given (t — t*) can vary 
by a factor of roughly 2. For example, a hundred thousand 
years after the first star forms, in one run the total stellar 
mass is WMq while in another it is 15M0. 


Cl Radial and Lateral Components of the 
Random Motion Velocity 

Figure C2 shows the radial vt, r and lateral vt, i = (vt, e + 
vt,4 >)C ^ components of the random motion velocity for 
the same collapsing region as shown in Figure (5), where 
Vt, e and vt, <j> are the random motion velocities along the 
9 and (f) directions defined from the z-axis. In the absence 
of self-gravity, a turbulent hydrodynamic cascade to small 
scales tends towards equipartition, (ut, r ~ vt,i), with a 
scaling behavior vt ~ similar to that seen in Larson’s 
size-linewidth relation; this is what we see in non-collapsing 
regions in our simulation. 

Figure C2 shows that both the radial and transverse 
components of the random motion velocity decrease with 
decreasing r for 0.4 < r < 3pc (except for a spike at 
r « 0.6pc). Furthermore, the ratio VT,r/vT,i ~ 1. We in¬ 
terpret the decrease as the decay of turbulence down a cas¬ 
cade. However, the decrease in both the total and in the 
longitudinal component, when fit with a simple power law, 
gives Vt ~ with p = 0.2, while the decrease in the ra¬ 
dial component of the turbulence corresponds to p ~ 0.35. 
Since both exponents are less than the value p « 0.5 that 
we see on larger scales or away from collapsing regions, we 
conclude that adiabatic heating is affecting both the radial 
and transverse components of the turbulence. 


^ Note that we define vt, i as an average so that we can compare 
it directly to vt, r- 
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Figure Bl. Histogram of the cumulative normalized M (left) and cumulative normalized solid angle (right) as a function of normalized 
density for r = 0.05 and 0.5 pc when the star reaches 1 and 4 Mq. 



Figure B2. Histogram of the cumulative normalized averaged M 
as a function of the cumulative normalized averaged solid angle 
for r = 0.05 and 0.5pc when the star reaches 1 and 4 Mq. 

At smaller radii, 0.04 pc r < 0.4 pc, the inward de¬ 
crease of both vt, r and vt. i slows and then reverses, as the 
flow passes r,. However, the ratio vt, t/vt, i is now only 1/2. 
Finally, at and inside the disk radius rd ~ 0.02 pc, the lat¬ 
eral turbulence once again decreases inward, while the radial 
component grows until much smaller radii, before decreasing 
again. 

If adiabatic heating is responsible both for the slower 
than normal decrease with random motion velocity at r > 
r*, and for the increase in random motion velocity inside r,, 
why does the ratio of the radial and lateral components of 
the turbulence vary? 

In Figure 5, \ur\ is decreasing with decreasing radius 
over the range 0.4 pc < r < 3 pc. What this decrease means 
physically is that as the gas falls in towards the center, it is 
being compressed not just in the 0 and (f) directions, but also 
radially. This compression along the radial direction should 
drive radial turbulence, while the lateral compression should 
drive lateral turbulence. This is why the radial and lateral 



Figure Cl. The run of velocity for a particle formed in a simula¬ 
tion using flash’s built in particle formation checks (as well as 
the density condition proposed by Padoan &; Nordlund (2011)). 
This shows that the dynamics of collapse are not overly sensitive 
to the star particle creation algorithm. 


components of the random motion velocity have the same 
magnitude. 

This physical reasoning also tells us that as the infall 
velocity increases inward over the range 0.04 pc < r < 0.4 pc 
(see Figure 5), the gas dilates in the radial direction even 
as it continues to compress in the transverse {9 and <j)) di¬ 
rections. Compression in the 9 and (j) directions will tend to 
drive an increase in the lateral components of the random 
motion velocity, but dilation in the radial direction will tend 
to drive a decrease in the radial component; of course both 
tendencies have to compete with (or add to, in the case of 
radial motion) the usual tendency for turbulence to decay, 
and the tendency, mentioned above, for hydrodynamic tur¬ 
bulence to tend to equipartition as the motion cascades to 
small scales. 

We interpret the rapid inward decline of VT,r starting 
at r « 0.4 pc as the effect of adiabatic cooling. The result is 
that the ratio vt.tIvt.i « 1/2 for 0.04 pc < r < 0.4 pc. 

Between « 0.02 pc and the local maximum of |Mr| 
at r « 0.04 pc, the infall velocity is large but roughly con- 
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Figure C2. The radial tix.r (blue solid), lateral ux,! = (ux,e + 
ux, c(l )/2 (red dashed) and total (green line over-plotted with dots) 
random motion velocities as a function of radius for particle B 
100, 000 yrs after the particle formed. The sound speed Cs is shown 
for comparison (black horizontal line). At large radii ux, r ~ ll 
for 0.04pc < r < 0.4pc ux,r < ux,i, while inside of 0.03pc ux,r 
quickly recovers and then exceeds ux i. The behavior of the lateral 
and radial velocities is dictated by the radial infall velocity in 
Figure 5 (see the main text). 

slant, meaning that the radial dilation ceases. We interpret 
the uptick in dt, r toward small radii as the result of the cas¬ 
cade of vt, i driving the radial component, as the turbulence 
strives to reach equipartition, combined with the cessation 
of adiabatic cooling associated with the cessation of radial 
dilation. 

Inside r « 0.02 pc the infall takes place primarily 
through a rotationally supported disk, in which both the 
vertical vt, e and azimuthal vt, </. component of the turbu¬ 
lence is greatly reduced (although we do not show the sepa¬ 
rate components in the figure). At the outer edge of this disk 
we see a sharp rise in the radial component of the random 
motion velocity, followed at yet smaller radii by a decrease 
in the total turbulent velocity. We interpret the drop in the 
total turbulent velocity at small radii as the flow settling 
into more ordered motion in an accretion disk. 


APPENDIX D: THE INITIAL MASS 
FUNCTION 

For completeness we report the IMF in this subsection, 
though we caution the reader again that, because we do 
not handle the thermal physics properly, the location of the 
break in the IMF is unlikely to be correct; however it is com¬ 
monly believed that the slope at the high mass end is set by 
the turbulence so that the thermal properties will not have 
much effect there. Figure D1 shows the IMF at the end of 
our Ramses run with Nj = 32. The plot includes a total of 
90 stellar particles, with a total mass of 24OM0, or about 
~ 1% of the total mass in the simulation. The time that 
Figure D1 is plotted corresponds roughly to the right edge 
(« 600, 000 years after the first star forms) in figure D2. It 
shows a form that is roughly consistent with observed IMFs, 


Figure Dl. Initial mass function for the Ramses run with Nj = 
32. The peak of the IMF is at 2Mq and the high mass slope is 
F = 1.36, where the Saltpeter slope is F = 2.35. The red line 
shows the least-squares fit to the mass function for M > I.OMq. 
Our IMF varies with time with the peak mass moving to lower 
mass and the value of F increasing with time. 



t-t* (yrs) 


Figure D2. Average stellar mass as a function of time for the 
Ramses run with Nj = 32 in green. Average stellar mass starts 
at Ri IMq increasing to r: 6 M 0 and then decreasing as low mass 
star formation begins. We also show the average stellar mass for 
Nj = 8 (blue) and Nj = 4 (red). This demonstrates that Nj = 4 
runs are converged. 

in that it has a power law at high masses, a peak around a 
solar mass, and a fall off at lower mass. The peak however, is 
at 2Mq which is about a factor of four higher than observed 
IMFs, and the fall off at high mass is too flat, indicating 
that we are top heavy. If the Saltpeter slope is denoted by 
F = 2.35, our slope is F = 1.36. 

Figure D2 shows the average stellar mass as a function 
of time. We see that the average mass is significantly higher 
than that of observed IMFs, where it is in the range of 0.3 — 
1.0. In addition, we see that this average mass rises initially 
as the massive stars grow and then decreases as low mass 
star formation kicks in. 
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Figure El. Plot of M*(i — t*) for Ramses runs with Nj = 32 
(blue), Nj = 8 (green), Nj = 4 (red). At the end of the Nj = 32 
run, the total stellar mass was M* Ri 24OM0. 

APPENDIX E: CONVERGENCE WITH Nj 

In this appendix, we examine how our results for the mass 
accretion rate depend on the resolution of the Jeans length, 
as quantified by Nj. 

Figure El shows the total mass in stars plotted as a 
function of time since the time t, at which the first star 
particle formed, for Nj = 4, 8, and 32. For the Nj = 32 
run, the final value of the total stellar mass M* « 24OM0, 
for Nj = 32) is about 1% of the total gas mass in the box, 
while the final t — t, ^ 0.6Myrs, about 15% of the free fall 
time 3.8Myrs for the mean density of the box. The green 
line shows the total stellar mass for Nj = 8, while the red 
line shows the same quantity for Nj = 4. 

The figure shows that for t — > 200,000 years the 

stellar mass as a function of t — t* is converged to within 
10 %, and to even better accuracy at late times. 

We have also done a convergence study for the average 
mass, see figure D2, showing that the mean stellar mass is 
converged for Nj = 4. This is consistent with the IMF being 
converged, albeit to a form that is not in good agreement 
with observations. We remind the reader that because of our 
use of an isothermal equation of state, we do not expect the 
IMF to match measured IMFs. 

As a further convergence check. Figure E2 shows the run 
of density as a function of radius for three different Ramses 
simulations. The Nj = 32 (blue) run had 3 star particles at 
0.5Mq < M„ < 3M0, while both the Nj = 8 (green line) 
run and the Nj = 4 (red line) had 9 star particles. We see 
convergence for all radii larger than the disk radius, rj- This 
illustrates that the density approaches an attractor solution 
that is robust against the underlying numerical technique. 

This paper has been typeset from a TEX/IAIl^X file prepared by 
the author. 


Figure E2. Plot of p(r) for Ramses runs with Nj = 32 (blue), 
Nj = 8 (green), Nj = 4 (red). These are the averaged profiles 
of the density for star particles with O.5M0 < M, < 3Mq. Each 
simulation had ss 88Mq worth of gas in star particles. The density 
for r > is the same all three runs, showing that the Nj = 4 
run is converged for r > rj. 
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